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ABSTRACT 

Feedback from massive stars is one of the least understood aspects of galaxy formation. We perform 
a suite of vertically stratified local interstellar medium (ISM) models in which supernova rates and 
vertical gas column densities are systematically varied based on the Schmidt-Kennicutt law. Our 
simulations have a sufficiently high spatial resolution (1.95 pc) to follow the hydrodynamic interactions 
among multiple supernovae that structure the interstellar medium. At a given supernova rate, we find 
that the mean mass- weighted sound speed and velocity dispersion decrease as the inverse square root 
of gas density. The sum of thermal and turbulent pressures is nearly constant in the midplane, so 
the effective equation of state is isobaric. In contrast, across our four models having supernova rates 
that range from one to 512 times the Galactic supernova rate, the mass-weighted velocity dispersion 
remains in the range 4-6 km s^^. Hence, gas averaged over ~100 pc regions follows P cx with 
a « 1, indicating that the effective equation of state on this scale is close to isothermal. Simulated 
H I emission lines have widths of 10-18 km s~^, comparable to observed values. In our highest 
supernova rate model, superbubble blow-outs occur, and the turbulent pressure on large scales is >4 
times higher than the thermal pressure. We find a tight correlation between the thermal and turbulent 
pressures averaged over '^100 pc regions in the midplane of each model, as well as across the four 
ISM models. We construct a subgrid model for turbulent pressure based on analytic arguments and 
explicitly calibrate it against our stratified ISM simulations. The subgrid model provides a simple yet 
physically motivated way to include supernova feedback in cosmological simulations. 
Subject headings: galaxies: formation — ISM: kinematics and dynamics — ISM: structure — methods: 
numerical — turbulence 



> 
m 



oo 

o 



X 



1. INTRODUCTION 

The observed vertical distribution of gas in the Galaxy, 
especially the cold neutral H I component at large 
heights, cannot be supported by thermal pressure alone 
(Lockman & Gehman 1991; Ferrara 1993). Often invoked 
to provide the required additional pressure are turbulent 
motions of the interstellar gas, magnetic fields, and cos- 
mic rays (Boulares & Cox 1990; Lockman & Gehman 
1991). Of these, turbulent pressure is perhaps the dom- 
inant pressure component in the interstellar medium 
(ISM; McKee 1990). There are strong reasons to believe 
that supersonic turbulence also regulates star formation 
processes in molecular clouds by lowering the efficiency 
of star formation (Elmegreen & Scalo 2004; Mac Low 
& Klessen 2004; Krumholz & McKee 2005; Krumholz & 
Tan 2007). 

Supernova (SN) explosions are thought to be the dom- 
inant agent for producing turbulent motions in star- 
forming parts of galaxies (Norman & Ferrara 1996; Mac 
Low & Klessen 2004). Interactions between blast waves 
from multiple SNe lead to compression and shearing of 
the gas. In extreme cases, explosions that are spatially 
and temporally correlated may collectively generate a su- 
perwind, observed at both low (e.g., Martin 1999; Heck- 
man et al. 2000) and high rcdshifts (Pettini et al. 2001, 
2002; Shapley et al. 2003; Weiner et al. 2008). McKee & 
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Ostriker (1977) provided a picture of the ISM in which 
three distinct phases regulated by SN explosions exist 
in rough pressure equilibrium. The multiphase model 
does not directly include turbulent pressure produced by 
interactions among multiple SN-driven blast waves and 
large pressure fluctuations in the ISM. Numerical simula- 
tions (Rosen & Bregman 1995; Korpi et al. 1999; Gazol- 
Patiiio & Passot 1999; de Avillez & Berry 2001; Kim 
2004; de Avillez & Breitschwerdt 2004, 2005; Mac Low 
et al. 2005; Slyz et al. 2005; Joung & Mac Low 2006, 
hereafter Paper I; Tasker & Bryan 2006) and observa- 
tions of thermally unstable gas (Hcilcs 2001; Hciles & 
Troland 2003, 2005) also support this dynamic picture. 
Since turbulent pressure is a significant, sometimes the 
dominant, form of pressure in the ISM, its inclusion is 
vital for any realistic model of the ISM. 

Another motivation for this study comes from long- 
standing problems in cosmology. The standard cold dark 
matter paradigm has been widely successful in explain- 
ing structures on cosmological scales (Cen et al. 1994; 
Zhang et al. 1995; Ranch 1998; Bertschinger 1998; Cen & 
Ostriker 2000; Springel et al. 2005; Dunkley et al. 2008). 
Yet our understanding of the luminous matter on galactic 
scales remains seriously incomplete, mainly due to com- 
plexities of baryon physics, for which additional nonlin- 
ear thermal and hydrodynamic processes must be taken 
into account. A prime example is the galaxy luminosity 
function, which shows lower abundances of galaxies both 
on the faint and the bright ends compared to the halo 
mass distribution (e.g., Benson et al. 2003). Other out- 
standing astrophysical problems include incorrect predic- 
tions for the size of disk galaxies (the "angular momem- 
tum problem"; e.g., Navarro & White 1994; Steinmetz & 
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Navarro 1999) and for the fraction of gas that should have 
cooled to form galaxies (the "overcooling problem"). 

Feedback from massive stars and active galactic nu- 
clei plays a central role to explain these discrepancies in 
recent semi-analytic models (Croton et al. 2006; Bower 
et al. 2006; Hopkins et al. 2008; Somerville et al. 2008). 
Despite mounting evidence that SNc played a major role 
in galaxy formation (sec, e.g., Adelberger et al. 2003 and 
references therein), it has been difficult to incorporate 
the details of stellar feedback into cosmological simula- 
tions. Early studies found that thermal energy input did 
not affect the ISM significantly, as it was radiated away 
almost instantly (Katz 1992; Navarro & White 1993). 
This was due to a lack of numerical resolution, which led 
to an absence of hot gas that would have radiated inef- 
ficiently and substantially decreased the overall cooling 
rate. Past attempts to resolve this difficulty (e.g., Kauff- 
man et al. 1993; Somerville & Primack 1999; Efstathiou 
2000; Thacker & Couchman 2001; Springel k Hernquist 
2003) have been phenomenological in nature, as they con- 
tain one or more free parameters that need to be adjusted 
to match observations, e.g., the observed rates of star 
formation. What is lacking is a detailed high-resolution 
3D model of turbulent flows within a small region of the 
ISM, to help identify basic properties of interstellar tur- 
bulence. 

Our aim in this paper is to provide a framework in 
which to handle SN feedback in cosmological simula- 
tions. We divide this challenging task into two parts. 
First, we run a suite of numerical simulations of strati- 
fied ISM driven by discrete SN explosions. The SN rate 
and the vertical gas column density are systematically 
varied over wide ranges to examine how the velocity dis- 
persion changes as a function of the gas density, length 
scale, and the assumed SN rate. In particular, we study 
the distribution and variation of turbulent pressure in 
these models. These local, high-resolution simulations 
have high enough spatial resolution to follow turbulent 
motions of the; interstellar gas in detail. 

The second step is to use the local ISM models to find a 
prescription for the physics unresolved in global, cosmo- 
logical simulations so that we have a better idea of how 
to handle SN feedback. We propose a subgrid model 
based on heuristic analytic arguments and simulation 
results. The idea is as follows: If we average relevant 
physical quantities over scales greater than the largest 
energy containing scale, we may describe the large-scale 
dynamics of the ISM that is insensitive to the details 
of the small-scale physics. This is motivated by similar 
studies in the past (Yepes et al. 1997; Springel 2000) as 
well as the new finding in Paper I that >90% of the total 
kinetic energy of the disk is contained on scales below 
200 pc, comparable to the size of resolution elements in 
current cosmological simulations (e.g., Governato et al. 
2004, 2007; Cen et al. 2005; Naab et al. 2007; Gibson 
et al. 2008; Joung et al. 2008). This suggests the in- 
teresting possibility that the SN-driven turbulence may 
be handled separately as subgrid physics in those simu- 
lations. Based on the physical characteristics of our lo- 
cal ISM models with discrete SN explosions, the subgrid 
model is physically well-motivated. We apply the results 
of our ISM models by formulating an accurate feedback 
prescription after deriving an effective equation of state 
for turbulent pressure of the medium. Such a prescrip- 



tion may be used to examine the effect of SN feedback 
in large cosmological volumes. 

Dib et al. (2006) investigated the relationship between 
the velocity dispersion of the gas and the SN rate in 3D 
simulations of unstratified, periodic boxes with a fixed 
mean density. Their aim was to explore the observed 
constancy of the velocity dispersion in the outer parts 
of galactic disks (Dickey et al. 1990; Kamphuis 1993) 
and the transition to the starburst regime. They sim- 
ulated a (1 kpc)^ volume of the unstratified ISM with 
128^ cells for various SN rates, while keeping the gas 
density of the box constant. They found nearly constant 
velocity dispersions for the H I gas at relatively low SN 
rates (0.01 to 0.5 times the Galactic value; note thtat 
they took a Galactic SN rate per volume of 2.58 x 10^ 
Myr~^ kpc""^), as well as a transition to the starburst 
regime at about the star formation rate per unit area of 
0.5-1x10"^ M0 kpc~^ yr~^ (assuming a Salpeter initial 
mass function; Salpeter 1955), in good agreement with 
observations. The present study extends their work to 
stratified atmospheres at higher spatial resolution where 
the gas density is systematically increased with increas- 
ing SN rate. 

The basic features of our ISM models are described 
in § 2. In § 3, we describe the thermal and structural 
properties of the models. We demonstrate in § 4 that 
the turbulent pressure of a medium with a given SN rate 
is independent of the local density, and study how the 
velocity dispersion and H I linewidth vary as a function 
of the assumed SN rate. Using these results, we construct 
a new subgrid model and calibrate it against the ISM 
simulations in § 5. Finally, we discuss several outstanding 
open issues (§6) and summarize our main results (§7). 

2. LOCAL ISM SIMULATIONS 

We use the model described in Paper I and its ex- 
tensions described here. Our simulations are performed 
using Flash v2.4, an Eulerian astrophysical hydrody- 
namics code with adaptive mesli refinement (AMR) ca- 
pability, developed by the Flash Center at the Uni- 
versity of Chicago (Fryxell et al. 2000). It solves the 
Euler equations using the piecewise-parabolic method 
(Colella & Woodward 1984) to handle compressible flows 
with shocks. For parallelization, the Message-Passing 
Interface library is used; the AMR is handled by the 
PARAMESH library. 

Although, in the present work, we consider SNe as the 
only source of interstellar turbulence, other processes 
such as gravitational instability (e.g., Li et al. 2005a), 
magneto- Jeans instability (Kim & Ostriker 2002), mag- 
netorotational instability (e.g., Piontek & Ostriker 2005), 
expanding H II regions (Matzner 2002) and protostel- 
lar outflows (Haverkorn et al. 2004; Nakamura & Li 
2007; Banerjee et al. 2007; Carrofl et al. 2008) should all 
contribute to the observed turbulence on various length 
scales. However, Mac Low & Klessen (2004) argue that 
these other processes may play relatively minor roles. We 
assume that they are insignificant on the spatial scales 
resolved in our models. 

Our computation box contains a volume of (0.5 
kpc)^x(10 kpc), elongated in the vertical direction. If 
the grid were fully resolved at our maximum resolution 
of 1.95 pc, it would contain 256^x5120 zones, but fully 
refined cells are mostly concentrated near the galactic 
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TABLE 1 

Parameters for the ISM simulations 



model 


S*,ii " 


^gas 


r " 




{risn) 


Ix 


11.0, 16.4 


1.87(6) 


8.50(-26) 


40 


14.8 


8x 


98.0, 147. 


8.23(6) 


1.55(-25) 


320 


27.4 


64x 


794., 1190. 


3.63(7) 


2.78(-25) 


2560 


41.4 


512x 


6360., 9540. 


1.61(8) 


5.11(-25) 


20480 


55.9 



Note. — The SN rates per unit area are in units of Myr 
kpc~^. Exponents are given parenthetically. 

^ Rates for (field, association) Type II SNe 

The gas column density at t = is given in Mq kpc"'^ 
The diffuse heating rate is expressed in erg 



midplane (the central ±200 pc) of our models. The 

gravitational potential of Kuijken & Gilmore (1989) is 
employed. For simplicity, the gravitational potential is 
fixed for all models, except in the model with the high- 
est SN rate (512x), where the increased gas surface den- 
sity dominates the gravitational potential, as it exceeds 
the assumed stellar mass surface density of 4.46 x 10'' 
M0 kpc"2 (Binney & Tr cmaine 1987). Assuming the 
Schmidt-Kennicutt law and taking Stot = 75 Mq pc^^ 
and Sgas = 5.3 Mq pc~^ for the solar neighborhood (Bin- 
ney & Tremaine 1987), we increase the gravitational ac- 
celeration due to the disk component in the Kuijken & 
Gilmore potential by a factor of 11 for the 512x model. 
This is our standard 512x model, and the results are 
based on this model, unless otherwise noted. For com- 
parison, we also ran a model with the same SN rate but 
with the gravitational potential used in the lower SN rate 
models; this latter model will be denoted 512xL. Self- 
gravity of gas is not included. Outflow boundary condi- 
tions are used in the upper and lower surfaces parallel to 
the Galactic plane, while periodic boundary conditions 
are employed elsewhere. 

For each SN explosion, we dump thermal energy Esn = 
]^q5i Qj-gg a small sphere whose radius varies as a func- 
tion of the local density. The radii are chosen such that 
radiative losses inside the spheres are negligible in the 
first few timesteps after explosions (Paper I); they usu- 
ally range between ~7 and ~50 pc. We use them to 
trace the thermal history of metal particles and to esti- 
mate their escape fraction from the galaxy. We assume 
that a fixed fraction, 3/5, of the Type II SNe are closely 
correlated in space and in time as a way of simulating 
super bubbles. The remaining explosions have random 
positions scattered through the galaxy mass, with scale 
heights of 325 pc for Type I and 90 pc for Type II SNe, 
to represent field SNe. The minimum number of SNe 
?^sn,min = 7, whilc the average and maximum numbers of 
SNe per superbubble (ngn) and ngn.max are determined 
by a power-law distribution (INb oc dn^ yielding the 
values given in Table 1 (Kennicutt et al. 1988; see also 
Paper I for details). We fix the Type I SN rate at the 
Galactic value = 6.58 Myr~^ kpc~^ in all models, 
as they take >1 Gyr to explode. 

We run a suite of numerical models with wide ranges 
of SN rates and gas column densities, as listed in Ta- 
ble 1. We adopt SN rates per unit area, S*, of 1, 8, 
64, and 512 times the Galactic rate (see Table 1 for the 
exact values adopted), and vary the gas column densi- 
ties systematically following the Schmidt-Kennicutt law 



(Kennicutt 1998). Wong & Bhtz (2002) reported ob- 
servational evidence for seven CO-bright spiral galaxies 
that the azimuthally averaged star formation rate per 
unit area scales with the gas column density similarly to 
the disk-averaged star formation law (see Li et al. 2006 
for similar results from numerical simulations). More re- 
cently, Kennicutt et al. (2007) reported evidence that the 
same star formation law is satisfied on scales of 0.5-2 kpc. 
The initial gas column densities Sgas are determined by 
inverting the relation 

where the coeffcient A — 2.5 x 10~^ Mq yr~^ kpc~^. As 
the Schmidt-Kennicutt law holds in various physical con- 
ditions found in galaxies, we believe that our choices of 
(Sgas, S*) are representative of typical regions in actual 
galaxies. 

Radiative cooling appropriate for an optically thin 
plasma with Z=Zq (Dalgarno & McCray 1972; Suther- 
land & Dopita 1993; see Fig. 1 of Paper I) and a diffuse 
heating term that accounts for the photoelectric heating 
(Wolfire et al. 1995, 2003) of the neutral gas are also in- 
cluded. The diffuse heating rate F in each model may not 
be directly proportional to because the far-ultraviolet 
(FUV) radiation responsible for the photoelectric heating 
{hv = 6-13 eV) may be significantly shielded in regions 
of high densities. In reality, the FUV radiation field is 
spatially nonuniform and time- varying (Parravano et al. 
2003). However, as radiative transfer is not included in 
our calculations, we resort to a simple scaling relation. 
To compute the mean intensity of the FUV radiation ap- 
propriate for each ISM model, we imagine a medium that 
emits and absorbs radiation uniformly. In the optically 
thick limit (appropriate for our box size of 500 pc), the 
intensity of the FUV radiation, /puv = jvloiv , where 
jn is the emissivity and = na is the absorption co- 
efficient with the number density of absorbers, n, and 
their cross-section, a. Assuming j,y oc oc Sg^g and 
n oc Egas (since the gas scale height is only weakly 
dependent on the SN rate, as shown later), we obtain 

/puv c< Egaf oc E*^*". With the additional assumption 
that F is proportional to the intensity of the FUV radi- 
ation, we obtain F cx E*' . As in Paper I, the diffuse 
heating rate in each model declines exponentially away 
from the galactic midplane with a scale height of 300 pc. 

We note that another set of models was run with 
a different scaling relation for the diffuse heating rate: 
r oc S* oc Sga^ . These models, however, yielded results 
inconsistent with observations. For example, almost all 
the gas in the 512x model in this set resided in the warm 
phase (T ~ 10''' K, the maximum gas temperature where 
diffuse heating is applied in our models), and the turbu- 
lent pressure was smaller than the thermal pressure by a 
factor of a few, in disagreement with observations. Our 
result is consistent with the simulations in Kim (2004), 
who found that the velocity dispersion of the gas in a 
SN-driven medium is lower for higher external pressure. 
Although we do not present these models in this paper, 
they illustrate the importance of including a proper form 
of diffuse heating rate, as it affects the balance among 
various interstellar gas phases. 
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3. THERMAL AND STRUCTURAL PROPERTIES 

Figure 1 displays (from top to bottom) gas density, 

temperature, and thermal pressure in the midplane of 
the four ISM models. Figure 2 shows probablity distri- 
bution functions for those variables. As in the Galactic 
SN rate model discussed in Paper I, the three classical 
phases of the ISM are present in every model, with cold, 
dense filaments surrounded by warm gases that in turn 
are embedded in hot media. However, their volume frac- 
tions vary as a function of the SN rate: The hot and 
cold phases increase in volume with increasing SN rate, 
while the volume occupied by the warm phase generally 
decreases (Fig. 2b) . Interestingly, the gas fraction in the 
thermally unstable range of temperatures (200 ^ T < 
7000 K) increases steadily with the SN rate (Fig. 2b), 
reflecting more dynamic ISM with higher SN rates. More 
subtle differences are also found. For example, the mean 
density of the hot gas becomes progressively higher as 
the SN rate increases (Fig. 2a), likely caused by an in- 
creased level of turbulent diffusion off dense shells and 
filaments where SN explosions occur more frequently (de 
Avillez & Mac Low 2002). Partly due to the increase 
in density, the overall thermal pressure Pther in the mid- 
plane increases with the SN rate (Fig. 2c). Within each 
ISM model, however, we find that thermal pressure varies 
within a range of roughly two orders of magnitude (Mac 
Low et al. 2005). The peak of the pressure distribution 
in the Galactic SN rate model lies at ~3 x 10'^ cm~'^ K, 
close to the mean value measured with C I fine-structure 
lines (Jenkins & Tripp 2007). Recent SN remnants arc 
associated with thermal pressures higher than the mean, 
whereas old SN remnants are tied to those that are lower. 
The dispersion in Pther rises for high SN rate models, 
where its maximum value reaches ~10'^ cm~^ K, close to 
the value measured near the center of M82 (Strickland 
& Heckman 2007). Our results are broadly in agreement 
with the analytic model by Monaco (2004) in which the 
ISM was modeled as a two-phase medium in pressure 
equilibrium. 

The vertical distributions of the ISM models are shown 
in Figure 3. In all models, correlated SN explosions 
create superbubblcs that vent hot (T > 10^ K) metal- 
enriched gas out of the galactic disk into the halo with 
high velocities {\vz\ > 300 km s~^). As the SN rate in- 
creases, the total gas mass also increases; a gradually 
larger mass fraction resides in the cold phase rather than 
in the warm phase. With stronger gravity, we would ex- 
pect the dense gas to be more confined to the galactic 
plane, decreasing the disk thickness, while superbubbles 
would blow out more easily into the halo. On the other 
hand, the change in the ISM near the galactic midplane 
should bo modest. This is what we find by comparing the 
two models, 512x and 512xL. The 512x model has higher 
gas densities in the midplane as well as higher turbulent 
velocities than those in the 512xL model. In particular, 
it is worth noting that the 512x model shows superbubble 
blow-outs, contrary to expectations from some previous 
work (see § 6 for more details). 

4. TURBULENT PRESSURE DISTRIBUTION 

We examine the distribution of turbulent pressures 

in our models after the systems reach statistical steady 
states at i « 80 Myr. Because kinetic energy is dis- 
tributed over a broad range of wavelengths (see Fig. 



8b of Paper I), turbulent pressure is inherently a scale- 
dependent quantity. To quantify the random (i.e., tur- 
bulent) component of the velocity field apart from the 
bulk motion of the medium, the root-mean-square (rms) 
velocity dispersion a is measured in the local frame of ref- 
erence, i.e., the center of mass frame of a selected volume. 
We choose the entire plane with \z\ < 125 pc, and tile it 
with small cubical boxes, starting from those with only 2 
zones (3.91 pc) on a side. The velocity dispersion within 
each box is computed from the random component of the 
velocity field w after subtracting the center-of-mass ve- 
locity of the box Vq from the original velocity field v. The 
box size is successively increased twofold until it reaches 
128 zones (250 pc) on a side. This procedure enables us 
to study the velocity dispersion as a function of scale. 
For each box of a given size, the turbulent rms velocity 
dispersion CTturb.b is calculated by 

^turb,b-Tj-/ PH'dV, (2) 

where Mb denotes the total mass within the box and 
w = V — Vq. Similarly, the thermal velocity dispersion 
(Tther,b is Computed by replacing \w\ with the local sound 
speed Cg. We prefer to use quantities directly related to 
pressure, so we set 

<7ther = (aLr,b/7)'/' , (3) 

aturb = (<7Lb,b/3)'/', (4) 

and 

f^tot = (<7ther,b/7 + 0-turb,b/3)^/^ (5) 

for ID thermal, turbulent, and total velocity dispersions, 
respectively (sec § 6). We hereafter drop the subscript 
"b" with the understanding that all velocity dispersions 
are scale-dependent quantities, computed for boxes of a 
particular size. Similarly, we define the total pressure 
Ptot as the sum of thermal (Pther = nkT) and turbulent 
(Pturb = /"^turb) prcssurcs. Bclow wc rcport our findings 
on the distributions of turbulent pressure found from this 
analysis. 

4.1. Pressure Equilibrium at Fixed SN Rate 

In Figure 4, thermal, turbulent, and total velocity dis- 
persions (Eq. 3-5) are plotted against average densities 
of individual boxes as a function of the box size. The 
plots on the right-hand-side display linear regressions of 
the corresponding plots in the left column, row by row. 
Small yellow circles represent the smallest boxes (3.91 
pc on a side), while big purple ones represent the largest 
boxes (250 pc on a side). The results are presented only 
for the Ix and 512x Galactic SN rate models at t w 80 
Myr, but we verified that the conclusions in this section 
hold for all four models throughout their steady-state 
evolution. 

The most striking feature of these plots is that, for 
low and intermediate densities {n < 10 cm~^) the mean 
velocity dispersion, a oc p^^^^. showing isoharic be- 
havior. This is true at every scale for both thermal 
and total velocity dispersion, and also at large scales 
for the turbulent velocity dispersion. The idea of ther- 
mal pressure balance between phases is not at all new 
(Field et al. 1969). Indeed, McKee & Ostriker (1977) 's 
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Fig. 1. — (From top to bottom) gas density, temperature, and thermal pressure in the midplane of our four ISM models. Each column 
corresponds to a single model as indicated at the top of each column: (from left to right), 1, 8, 64, and 512 times the Galactic SN rate. 
The colorbars shown on the right apply to all models. 



three-phase ISM model was also based on the premise 
of "rough pressure balance" between phases, motivated 
in part by Chandrasekhar & Fermi (1953) and Spitzer 
(1956). What is new and remarkable is that the turbulent 
pressure is likewise in approximate pressure equilibrium 
(Fig. 4b) on large scales where Pturb ^ ^thcr- In particu- 
lar, at any given scale, CTtot oc indicating that the 
total pressure Ptot = P <^tot nearly constant in the mid- 
plane. This implies that different gas tracers, which often 
probe narrow ranges of gas densities, may find discrepant 
kinematics (in particular, velocity dispersions) but that 
there may be an underlying correlation between them. 
The pressures are in approximate equilibrium, although 
the gas density and temperature individually vary by ^7 
orders of magnitude. 

Furthermore, the thermal velocity dispersion is scale- 
independent, as indicated by the lines in the top- left 
panel of Figure 4, which lie almost on top of one an- 
other. In contrast, the turbulent velocity dispersion in- 
creases with scale. As the box size increases, larger eddies 
gradually contribute to the turbulent velocity dispersion, 
and CTturb increases just as the power spectrum of kinetic 
energy predicts (§ 4.4). As a result, the thermal pres- 
sure dominates the total pressure on small scales, while 
the turbulent pressure dominates on large scales. At any 
particular scale, the scatter is about an order of magni- 
tude in (T, leading to ~2 orders of magnitude scatter in 
the pressure (see Fig. 7c of Paper I). 

4.2. Turbulent Velocity Dispersion vs. SN Rate 
4.2.1. Mass-weighted velocity dispersion 



How does turbulent pressure scale with the SN rate? 
To answer this question, we compare crturb of the sec- 
ond largest boxes with 125 pc on a side in all models 
(purple circles in Fig. 4b) and obtain relevant physical 
quantities averaged over 125 pc scales. (Hence, unless 
otherwise noted, trturb and Pturb henceforth refer to the 
values computed in boxes with 125 pc sides.) Using these 
models, we can find how turbulent pressure depends on 
the SN rate, or after conversion via an appropriate IMF, 
the star formation rate. 

Figure 5 displays the ID mass-weighted turbulent ve- 
locity dispersion against the SN rate for the four SN rates 
and gas column densities employed in our models. The 
mean velocity dispersions Uturb (in the mass-weighted 
rms sense) and their standard deviations are given in Ta- 
ble 2 The changes in the mean values are minimal across 
the range of 512 in SN rates: the best-fit line has the 
form 

atu.b = (5.6 ± 0.8) i]^-0.045± 0.033 ^-i (g) 

The variation of the mean total (i.e., thermal plus tur- 
bulent) ID velocity dispersions is shown in solid line, 
bounded by one sigma dispersions (dotted lines). For 
the mean values for this quantity given in Table 2 the 
3D velocity dispersions are 8-18 km s~^, comparable to 
the sound speed of the warm gas. The scatter in the 
velocity dispersion decreases with increasing SN rate. 

For a resolution study, we ran the 512xL model with 
spatial resolutions half (3.91 pc/cell; left) and double 
(0.98 pc/ccll; right) that of the fiducial model (1.95 
pc/cell; middle). The turbulent velocity dispersions mea- 
sured in these models are indicated by the error bars 
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Fig. 2. — (From top to bottom) Probability distribution func- 
tions of (a) gas density, (b) temperature, and (c) thermal pressure 
near the midplane {\z\ < 125 pc) of all models showing fractional 
volumes contained in logarithmic bins. The Ix, 8x, 64x, and 512x 
models are shown in solid, dotted, dashed, and dot-dashed curves, 
respectively. 



TABLE 2 

Results from the ISM simulations 



model 


0"turb 


Ctot 


CTv.H I " 




/w 


Ix 


6.1 ± 3.3 


10.6 ±5.9 


12.5 


137 


0.460 


8x 


4.7± 1.5 


7.3 ±2.4 


17.8 


87 


0.274 


64x 


4.0 ±0.5 


5.4 ±0.9 


12.4 


57 


0.431 


512x 


4.7 ±0.4 


5.2 ±0.5 


14.4 


55 


0.609 



Note. — See text for definitions of the variables. 
^ The first three columns are in units of km s~^. 
^ In pc 
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Fig. 3. — Vertical distributions of the four ISM models in this 
paper, (a) Ix, (b) 8x, (c) 64x, and (d) 512x, from a snapshot in 
time. For each model, (from left to right) slices of gas density, 
temperature, thermal pressure, and z-component of the velocity, 
and the projected gas column density are shown for the middle 2 
kpc of the models {\z\ < 1 kpc). 

in thick grey lines without open circles, offset along the 
X-axis for clear presentation. The velocity dispersions 
are 3.3 ± 0.6, 3.8 ± 0.4, and 4.0 ± 0.3 km s^^ for the 
low, intermediate and high resolution models, respec- 
tively. Therefore, although the 512xL models have not 
fully converged, the velocity dispersion is likely converg- 
ing to a finite value at infinite resolution, since factor of 
two changes in resolution led to only ~24% and ~5% in- 
creases in the velocity dispersion, with rapidly declining 
fractional changes. Since the 512xL model is associated 
with the highest mean density and hence most affected 
by resolution-dependent effects, we expect better conver- 
gence in the other three models. 

4.2.2. H I linewidths 

High-resolution H I emission maps of nearby, almost 
face-on galaxies are now available (e.g., Petric & Rupen 
2007; Tamburro et al. 2008). These observations show 
how the H I line widths and shapes change as a function 
of the galactocentric distance, gas density, distance from 
spiral structure, and star formation rate. Hence, they 
provide important constraints on the dominant driving 
mechanism of the neutral gas. 

The brightness temperature, T^, is computed by (e.g., 
Spitzer 1978) 



Tb (w) = (5.49 X 10"^^ K) 



H I 



A 



.w 



'1 



(7) 



where w is the velocity offset from the 21 cm line cen- 
ter, Nn I is the column density of H I in cm~^. Aw is 
the velocity width in km s^^ that we take to be 2.5 km 
s~^, which was the velocity resolution in Petric & Rupen 
(2007), and the velocity-dependent optical depth is given 

by 



T,„ = 5.49 X 10 



-19 



7ri/2 5T 



(8) 



where the Doppler lincwidth b — {2kT/m,)-^/^ in units of 
km s""'^ . Since we do not follow the multiple species of hy- 
drogen explicitly in our simulations, our estimate of A^h i 
is based on simple temperature cutoffs. Specifically, we 
assume that 100% of hydrogen between 50 and 7000 K 
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Fig. 4. — Mass-weighted (a) thermal, (b) turbulent, and (c) 
total velocity dispersions plotted against average densities of the 
boxes. The data are taken near the midplane of our Galactic SN 
rate model for various box sizes over which the velocity dispersions 
were computed. Small yellow corresponds to the smallest cubic 
boxes with 3.91 pc on a side, and purple to the largest boxes with 
125 pc on a side. The plots on the right-hand-side display linear 
regressions of the corresponding plots in the left column, row by 
row. On a log-log plot, the distribution of points can be character- 
ized by a slope of —1/2, implying a nearly constant total pressure. 
In (b) and (c), the velocity dispersion increases with scale because, 
as the box size increases, larger eddies gradually contribute to the 
turbulent velocity dispersion, (d)— (f) are similar plots for the 512x 
model. 



Fig. 5. — Turbulent ID velocity dispersions in boxes with 125 pc 
sides near the midplane (open circles), plotted against the SN rates 
in our four models. The changes in the mean values (represented 
by the X's and error bars) are minimal across the range of 512 in 
SN rates. The points in thick grey lines indicate turbulent velocity 
dispersions measured in the 512xL models with spatial resolutions 
lower (3.9 pc/cell; left) and higher (0.98 pc/cell; right) than the 
fiducial model (1.95 pc/cell; middle), offset along the x-axis for 
clear presentation. The variation of the mean total (thermal plus 
turbulent) ID velocity dispersions are shown in solid line, bounded 
by one sigma dispersions (dotted lines). The grey dashed lines 
represent the linewidths of H I emission, computed with the entire 
line profiles (upper curve) and with only the line centers (lower 
curve), as described in the text. 

exists as H I, while the neutral fraction increases linearly 
from 30 to 50 K. Below 30 K, hydrogen is assumed to be 
fully molecular. Note that the H I linewidths are gener- 
ally larger than the turbulent velocity dispersions. This 
happens because the bulk of the highest density g as , as- 
sociated with the lowest velocity dispersion, is molecular 
and so excluded in the calculation of the former quantity. 
The uncertainty in H II fraction and in H2 fraction affect 
the line profiles in the broad wings and in the centers, 
respectively. 

Figure 6 displays the 21 cm H I emission line profiles 
for the four ISM models as viewed face-on. With the 
exception of the 8x model, the line shapes are similar to 
the universal profile found by Petric & Rupen (2007). As 
shown, the profiles are wider than single Gaussian fits. 
The broad wings may be due to the warm component of 
the neutral gas. 

We estimate the linewidth, (Tv,h i, by fitting a single 
Gaussian to each line profile (grey dashed curves). The 
Gaussian fits poorly match the line profiles in all mod- 
els (although it is particularly bad for the 8x model), 
because of the broad wings, so we additionally fit only 
the central ±15 km of the line profiles with nar- 
rower Gaussian functions shown in grey dotted lines. 
The broad wings in the 8x model are due to warm 
gases that are more prominent at lower temperatures 
(4000 < T < 7000 K) than in the other models (see 
Fig. 2b). If we set the maximum cutoff temperature 
for neutral gas to 4000 K instead of 7000 K, the wings 
show substantially decreased flux as in the other models. 
On the other hand, removing the gas at high altitudes 
{\z\ > 0.25 kpc) does not affect the line profile, which in- 
dicates that the broad wings are not caused by a galactic 
wind. 
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Fig. 6. — Hi line profiles in the four ISM models viewed from 
the top. {Top row) Ix, 8x, {bottom row) 64x, and 512x models. 
Single Gaussian function fits, shown in grey dashed curves, poorly 
match the line profiles, particularly for the 8x model, because of 
significant broad wings. Hence, we additionally fit only the central 
±15 km s~^ of the line profiles with narrower Gaussian functions 
shown in grey dotted line. The H I linewidths measured near the 
line centers change relatively little as the SN rate goes up. The 
tiny dip near the emission peak in the 512x model is caused by 
self-absorption of H I emission. 

We find that the Unewidths (fit to the central ±15 km 
s~-^) crv,H I are nearly constant across a wide range of SN 
rates as given in Table 2, similar to the behavior of the 
mass-weighted velocity dispersions. These linewidths are 
comparable to the velocity dispersions found by Agertz 
et al. (2008) in global galaxy models. They attrilnited 
the turbulent motions to global non-axisymmctric modes 
and cloud-cloud tidal interactions and merging. We find 
that SN explosions drive turbulent motions that are at 
least comparable in terms of energy. The H I linewidths, 
associated with only neutral hydrogen, are systemati- 
cally higher than the mass- weighted velocity dispersions. 
Hence, using H I linewidths to infer the overall velocity 
dispersion may lead to an overestimate. Our results are 
consistent with the observed near constancy of line-of- 
sight velocity dispersions with the galactoccntric radius 
(Dickey et al. 1990; Kamphuis 1993), and extend the 
numerical results of Dib et al. (2006) to stratified inter- 
stellar media with gas surface densities increasing with 
SN rate. 

4.3. Thermal vs. Turbulent Pressures 

Figure 7a shows a remarkably tight correlation between 

thermal pressure and turbulent pressure (computed in 
boxes with 125 pc sides) across our ISM models. This 
arises because higher SN rates naturally lead to increases 
in both thermal and turbulent pressures. The best-fit 
line has a slope that is steeper than unity (1.97 ± 0.13), 
implying that the turbulent pressure increasingly domi- 
nates at higher SN rates. Note that this does not nec- 
essarily mean that turbulent pressure is comparable in 
magnitude to thermal pressure on all scales. Turbulence 
from large-scale motions (i.e., scales larger than our box 
size that are not well modeled by the simulation) may 
have a significant contribution, and so Ptot / -Pther niay be 
significantly larger than ~2 on scales above 125 pc. 



4.4. Kinetic Energy Spectrum 

In Paper I, it was shown that the kinetic energy in 

an explosion-dominated medium is distributed over a 
broad range of wavelengths. Then how do the charac- 
teristic scales of the density and kinetic energy distribu- 
tions change as a function of SN rate? Figure 8 displays 
angle-averaged density power spectra (left) as well as ki- 
netic energy spectra in black and angle-averaged veloc- 
ity power spectra in grey (right); see Paper I for details 
on how they are computed. As in the Galactic SN rate 
model, the velocity power spectra deviate significantly 
from the Kolmogorov spectrum, which is applicable to 
an incompressible medium driven at a single large scale: 
these assumptions are not applicable to the ISM driven 
by supersonic shocks on a range of scales. It is clear from 
Figure 8 that the scale containing most of the kinetic en- 
ergy shifts to larger wavenumbers (smaller wavelengths) 
with increasing SN rate. We quantify this change by 
computing the effective driving scale^ defined by 



L 



We find the values of id,eff given in Table 2. 
The best-fit line has a form id.eS ~ (126 ± 
18 pc) (SjS.^i) -"-isio.os^ ^Yiere S.^i is the SN rate 
per unit area in the Ix model. It should be remem- 
bered that we are dealing only with SN driven turbu- 
lence. Large-scale turbulence driven by, e.g., self-gravity 
may contribute to turbulence on larger scales, so that the 
combined power spectrum may be more or less continu- 
ous out to kiloparsec scales, as observed (e.g., Elmegreen 
et al. 2001). 

5. SUBGRID MODEL FOR TURBULENT PRESSURE 
5.1. Basic Features 

A genuine, three-dimensional, cosmological model that 
properly includes SN feedback would need to encom- 
pass an unrealistically large dynamic range, starting from 
length scales responsible for galaxy formation (>1 Mpc) 
down to those relevant to star formation (<1 pc). The 
resulting spatial dynamic range of at least 10^ requires 
computational resources currently impracticable. Past 
attempts to circumvent this problem include analytic 
approaches based on McKee & Ostriker (1977) 's mul- 
tiphase model of the ISM (e.g., Efstathiou 2000) and 
semi-analytic methods with parametrized prescriptions 
for star formation and feedback (Kauffmann et al. 1993; 
Somerville & Primack 1999; Cole et al. 2000). 

An alternative method is to use subgrid (or subresolu- 
tion) models (Ycpcs et al. 1997; Springel 2000; Semelin & 
Combes 2002; Springel & Hernquist 2003). In these mod- 
els, one represents the physics unresolved in the global 
model (i.e., turbulent motions on small scales) by an 
analytic prescription. Relevant physical quantities are 
evolved on a subgrid scale, and their appropriate aver- 
ages yield a simple description of the medium on the 
smoothing scale. The desired net effect is to have a large- 
scale simulation on a coarse grid, indistinguishable from 
the one in which all motions are resolved down to parsec 
scales. 

Note that the term "driving scale" is a misnomer since there 
is no single scale where energy is injected, unlike in Kolmogorov's 
idealized picture of incompressible turbulence. 
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Fig. 7. — Turbulent pressure variation in our four models with different SN rates. Displayed are: (a) turbulent vs. thermal pressure, (b) 
assumed SFR vs. density, (c) turbulent pressure vs. SFR, and (d) turbulent pressure vs. density. Each quantity is evaluated in (125 pc)"^ 
boxes near the midplane of the model galaxies. Different symbols indicate the ISM models; Ix (circles), 8x (triangles), 64x (diamonds), and 
512x (squares). In (a), the black dashed line is drawn for Pturb = -fther- The dark solid lines represent best-fit lines to all four ISM models. 
The grey dashed lines in (c) and (d) represent the relations expected from the subgrid model if the ratio ej/e2 is constant, equations (20) 
and (19), respectively. 

ergy of the gas, a term that helps to support the gaseous 
medium in the absenee of substantial thermal support. 
Nonlinear physical processes such as shock compressions, 
radiative cooling and thermal instability shape the prop- 
erties of interstellar turbulence. Hence, a subgrid model 
based on heuristic analytic arguments should be checked 
by direct numerical simulations that include the relevant 
physics. For this reason, we explicitly calibrate our sub- 
grid model against the results of our local ISM simula- 
tions. The subgrid model is therefore based on the results 
of our high-resolution local ISM models that resolve tur- 
bulent motions of the gas. In a sense, we imagine an 
entire model galaxy on a coarse grid (~100 pc spatial 
resolution), and within it, identify each cell with a box 
with side ~100 pc taken from one of our local models. 

The key parameter in our subgrid model is the turbu- 
lent pressure. Local ISM models described in this paper 
and elsewhere suggest that turbulent pressure is more 
important than either thermal (Boulares & Cox 1990; 
Lockman & Gehman 1991) or magnetic pressures (Kim 
2004; de Avillez & Breitschwerdt 2005) in most ISM con- 
ditions. This suggests that cosmological subgrid models 
should include a term for the turbulent pressure, rather 
than focusing on thermal pressure alone. Here we sug- 
gest a simple extension where we add an additional en- 
ergy term in the fluid equation to describe the turbulent 
energy on scales not resolved by the simulation (Joung 
2006; Bryan 2007). Two characteristics of our local ISM 
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Fig. 8. — Left: Angle-averaged density power spectra. Right: 
Kinetic energy spectra (black) and angle- averaged velocity power 
spectra (grey). The kinetic energy is distributed over a broad range 
of wavelengths, but the characteristic scale containing most of the 
energy shifts to larger wavenumbers (shorter wavelengths) as the 
SN rate increases. 



Our approach is inspired by the subgrid model de- 
scribed by Springel (2000) , which assumes that turbulent 
motions below resolved scales produced by SN feedback 
can be represented by a second reservoir of internal en- 
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models arc conducive to developing such a subgrid model 
for turbulent pressure as a way to represent small-scale 
motions. First, >90% of the total kinetic energy is on 
scales below 200 pc (Paper I). In current high-resolution 
cosmological simulations, therefore, most of the turbu- 
lent energy lies on scales smaller than a resolution ele- 
ment, and there is hope that such a model assuming a 
separation of scales would work. Second, the turbulent 
pressure is nearly constant on scales greater than ~100 
pc (§ 4.1), which suggests that for a given SN rate, a 
single value of turbulent pressure can characterize the 
medium regardless of the local gas density, albeit with 
the scatter shown in Figure 4b {left). 
In our subgrid model, the pressure is given by 

P = Pther + Pturb = (T " 1) P {<i + <?) (10) 

where e is the specific thermal energy density, while q = 
cr^mb is the specific turbulent energy density. Essentially, 
wc require an additional energy equation that is similar 
to the one for thermal energy but includes a source term 
reflecting the turbulent energy input from SNe instead 
of thermal heating terms, and a sink term corresponding 
to the decay of turbulence instead of radiative cooling 
terms: 



dq 
dt 



-(7 -l)qV -v + ei 



Efb P* 



£2- 



(11) 



P - L 

The star formation rate per unit volume, /9* = 
c*{p/tdyn), where tdyn = [37r/(32Gp)]^/^ is the dynami- 
cal time and the star formation efficiency per dynamical 
time, « 0.02 0.03 (Krumholz & Tan 2007). The local 
energy input due to SNe, ein = £tbP*, where the energy 
input per unit stellar mass, efb = 4 x 10^* erg Mq^ if the 
initial mass function (IMF) with a slope of —1.5 between 
0.1 and 40 M© and all stars with mass >8 Mq exploding 
as SNe. The turbulent energy is assumed to decay at a 
rate proportional to q/tcmss where Across = L/q^^^ is the 
crossing time at the turbulent driving scale (Mac Low 
1999); we assume L » id.eff) calculated in § 4.4. 

The two dimensionless parameters ei and €2 are se- 
lected to match local ISM simulations. The first of 
these parameters, ei, measures the efficiency of conver- 
sion of SN energy into kinetic motions, which may de- 
pend on, e.g. , the density and metallicity of the surround- 
ing medium (e.g., Thornton et al. 1998). We expect ei 
to decrease with increasing SN rate. The second param- 
eter, on the other hand, may be more or less constant; 
Mac Low (1999) found £2 — 0.84 for a range of energy 
injection rates and driving scales from a large set of hy- 
drodynamic and magnctohydrodynamic simulations with 
an isothermal equation of state. 

The thermal energy equation, in turn, is given by 

= -(7-l)eV-t;-h + (l-ei) +£2^-, 

dt p p L 

(12) 

where nT{p, e, Zisrf ) and n A{p, e) denote appropriate 
thermal heating and radiative cooling rates per unit vol- 
ume, respectively, n is the number density of the gas, 
and /isii is the intensity of the interstellar radiation field. 
The last two terms account for direct thermal heating 
from SNe and the dissipation of turbulent energy into 
heat. When combined with the Euler equations, equa- 
tions (10), (11) and (12) determine the evolution of ther- 
mal and turbulent energy densities. 



5.2. Calibration with Local ISM Models 

We can make simple analytic arguments to see how 

turbulent pressure is expected to vary as a function of the 
SN rate. To do this, we assume a steady state, where the 
energy input rate from SN explosions, Cin, is balanced by 
the dissipation of turbulent motions. We set dq/dt = 
and V • w = in equation (11) to obtain 

,3/2 



e'i(l - /w)em = £2 



pq 



(13) 



where we write ei as e[{l — /w) to explicitly account 
for the dependence on /w, the fraction of input en- 
ergy that goes into driving large-scale (> 100 pc) mo- 
tions such as galactic winds, and hence is not captured 
in the subgrid model. Now, to facilitate comparison 
with our ISM simulations, where we used the Schmidt- 
Kennicutt law for the initial conditions (Eq. 1), we adopt 
£in = efb^(2/5Hgas)^-''/(2HsN): where Fgas and Hsn are 
vertical scale heights for gas and SNe, respectively. The 
turbulent pressure is then given by 

-Pturb^pg = p 1 : — ) (14) 



£2 



p J 



c[(i-u)Hy!L]y'p''/'^ (15) 



where the combination of parameters, 

£2 -f^SN / 



(16) 



Now, consider each factor in brackets on the RHS of 
equation (14). The driving scale has already been com- 
puted (Eq. 9); we take L Ld.eff cx 



-0.15 ±0.03 



DC 



-0.15 ±0.03 



, if Efb and /fsN change little with the SN 

rate. To compute the gas scale height, -ffgas, we plot the 
SN rate (normalized to the Galactic rate) as a function 
of gas density in Figure 7b. If iJgas is constant across 
our models, the energy input rate will scale with the gas 
density in the same way as it varies with the gas surface 
density ( cx p^''^)- This is indeed close to what we 
find for the best-fit line shown in Figure 7b; 

\ 1.33 ±0.03 

P 



(1.79 ± 0.17) S.i 

\Pi 



(17) 



where pi ~ 0.7 cm~'^ is the mean density in the mid- 
plane of the Ix model. From this relation and ifsN ~ 90 
pc, we infer a nearly constant gas scale height of 170 pc 

(E,/E,,i)0-"4±o-02. 

To estimate how the fraction of energy driving large 
scale motions (/w) scales with the SN rate, we distin- 
guish between the small-scale, turbulent motions within 
boxes having 125 pc sides, as examined in § 4, and the 
bulk motion of these individual boxes. This amounts 
to comparing the sum of the kinetic energy contained 
in (mass-weighted) velocity dispersions inside the boxes 
(pcr^^j.,^/2) and that from the center-of-mass velocities of 
the boxes themselves {p\vq\'^/2). Using this procedure, 
we obtain the values of /^v given in Table 2. This does not 
contradict our earlier result that almost all the turbulent 
energy is contained on scales below ~200 pc, because we 
were then discussing the kinetic energy spectrum in the 
midplane region only. In fact, if we confine our atten- 
tion to the region with \z\ < 125 pc and compute the 
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energy fraction in large-scale motions, we obtain signifi- 
cantly smaller values (0.084, 0.065, 0.149, and 0.455). In 
other words, most of the kinetic energy on large scales 
comes from the bulk motions at high altitudes: galactic 
winds. If we fit a straight line through the values of 
given in Table 2, we find 



-0.06 ±0.05 



1-/^ = (0.65 ±0.15) 



(18) 



which is effectively a constant. 

To study the variation of turbulent pressure in our ISM 
simulations, we compute Pturb in boxes of size (125 pc)^ 
taken from the midplane region {\z\ < 125 pc) of the 
four ISM models. The results given in Figure 7c-d show 
the turbulent pressure depends on the input SN rate as 
well as the gas density averaged over cubes with 125 pc 
sides. The expectations from the subgrid model (Eq. 14) 
assuming constant C are shown as grey dashed lines in 
Figures 7c and 7d. The best-fit lines to the data, shown 
in dark solid lines, have the form 



turb 



oc E 



0.66 ±0.04 



(19) 

and, in terms of the gas density, 

Pturb OCpO•««±0•0^ (20) 

This effective equation of state (EOS) is slightly more 

compressible than an isothermal EOS, for which P p^. 
In contrast, much stiffer effective EOSs were adopted in 
previous subgrid models. For example, Springel (2000) 
used an EOS of the form P oc p^, obtained by assuming 
that all SN energy first goes into turbulent energy and 
then converts to thermal energy with an efficiency that 
depends specifically on gas density (ei = 1 and £2 oc 
p~^/'^ in our notation). In the hybrid multiphase model 
developed by Springel & Hcrnquist (2003), hot gas from 
SNe delays radiative cooling and provides the necessary 
additional pressure. Here, the effective EOS was also 
stiffer than isothermal above a threshold density ('-^O.l 
cm~^), where star formation was assumed to switch on. 
Our result indicates that SN-driven turbulence does not 
naturally provide a stiff EOS on 100-200 pc scales but 
rather that the effective EOS of the ISM averaged over 
these scales is slightly sub-isothermal (sec § 6 for further 
discussion). 

The simulation results deviate slightly from the naive 
expectations based on constant C, especially in the mod- 
els for starbursts (the 64x and 512x models). This sug- 
gests that C may be close to but not exactly a constant. 
Combining the scaling relations we have found, we obtain 

Coc£i-0•^5±0•°^ (21) 

or e'l oc ei;°-82±0-ii oc p-i-i5±o.i5 fo^ a constant £2- 
Hence, a smaller fraction of SN energy is used to drive 
kinetic motions as the SN rate (and the mean gas density) 
increases, although the dependence is weak. 

We must mention several caveats that may affect our 
conclusion. Our ISM simulations do not include mag- 
netic fields or self-gravity of gas, which can change the 
structure and the pressure distribution of the medium, 
especially in the dense regions. In principle, unphysical 
cooling of numerical origin may be important, particu- 
larly in the 64x and the 512x models. This would occur 



at the interfaces between hot and cold gas, which are in- 
evitably wider than the physical value, and so contain 
more gas at intermediate temperatures between 10^ K 
and 10^ K that is subject to strong radiative cooling than 
they physically should (e.g., Mac Low et al. 1989). Such 
cooling could yield values of Pturb that are lower than the 
correct values. However, the resolution study discussed 
in § 4.2 shows that this is a surprisingly minor problem 
in practice. An alternative method to examine this is- 
sue would be to use a tracer field (Yabe & Xiao 1993) 
to suppress the excessive cooling in yoimg SN remnants 
as in previous work (Mac Low & Ferrara 1999; Fujita et 
al. 2004), after appropriately modifying the criterion for 
suppression for the case of multiple discrete SN explo- 
sions. 

6. DISCUSSION 

Comparison with Chandrasekhar's formula for the effective 
sound speed. — In § 4, we have used Ptot = p(o'ther + 
o-Lb) = PCo-uicr.b/T + o-turb,b/3)- Thc reader may won- 
der why this is the appropriate expression. For example, 
compare it with the expression c^^g = 0"^ + (w^)/3 where 
Cg^off is the effective sound speed, which Chandrasekhar 
(1951) used to extend Jeans's analysis to an infinite ho- 
mogeneous turbulent medium. One immediately notices 
that, in his formulation, Cg was used instead of Cs/7. To 
understand this, we propose the following thought exper- 
iment. If we have two parcels of gas at the same temper- 
ature but with two different adiabatic indices (e.g., 7 = 
5/3 and 7/5), they slioidd have the same pressure (since 
Pther = nkT). However, when you try to compress the 
gases, they have different resistance to the compression, 
i.e., the one with higher 7 is more difficult to compress 
and hence their sound speeds differ. The first case is rel- 
evant to pressure equilibrium in our models; the second 
to Chandrasekhar's problem. Essentially, our problem 
is a pressure issue, while Chandrasekhar's is a commu- 
nication timing issue (the Jeans stability criterion says 

''Cross < ^dyn)- 

Large-scale turbulence driving. — We argued that most of 
the kinetic energy due to SN-driven turbulence can be 
contained in resolution elements (with size of 100-200 
pc) of cosmological simulations, hence allowing the ap- 
proach of subgrid modeling. In reality, however, there 
are numerous other energetic processes that occur on 
larger scales that may affect the dynamics. An exam- 
ple would be gravitational instability induced by spiral 
structure such as the magneto- Jeans instability (Kim & 
Ostriker 2002). Li et al. (2005a, 2006) argued, using an 
isothermal equation of state to represent stellar feedback, 
that gravity controls when and where star formation oc- 
curs, and reproduced the Schmidt-Kennicutt law with 
realistic star formation thresholds. The use of a subgrid 
model for interstellar turbulence does not mean that we 
ignore these larger-scale, possibly gravity- induced mo- 
tions. Rather, we are making an assumption that the 
characteristic lengths at which these two processes op- 
erate are sufficiently far apart that they can be treated 
separately in global models. One way to test this hy- 
pothesis is to compare simulated and observed column 
density power spectra (Elmegreen et al. 2001) for both 
local and global simulations of galaxies. 
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Observed linewidths in ULIRGs. - High velocity disper- 
sions in molecular gas, crco ~ 100 km s^^, were observed 
in ultraluminous infrared galaxies (ULIRGs) such as Arp 
220 (Downes & Solomon 1998). Murray et al. (2005; see 
also Thompson et al. 2005) suggested that in extreme 
starbursting systems such as ULIRGs, the ambient high 
density gas impedes the expansion of SN remnants, dras- 
tically lowering the efficiency with which SN energy es- 
capes the host galaxy. Their model does not include 
the effect of superbubbles: correlated SNe with high nsN 
may decrease the ambient density near explosion sites, 
and the overall radiative cooling rate in such explosions 
may be significantly reduced. In our models with high 
SN rates, SN explosions near the midplane drive large- 
scale motions about as effectively as in models with lower 
SN rates. The fraction of energy stored in turbulent mo- 
tions as opposed to heat increases with SN rate. We 
note that Fujita et al. (2008) recently found that ex- 
tremely superthermal linewidths (~320 km s~^) in Na I 
absorbing gas could be reproduced using a model for an 
energy-driven superbubble with a spatial resolution of 
< 0.2 pc, assuming a single central starburst region and 
axisymmetry. The energy input rate of our 512x model is 
~1.3xl0''^ erg s~^, comparable to the lowest mechanical 
luminosity model in Fujita et al. (2008). 

Effective equation of state. — The models presented in this 
paper demonstrate that a stiff EOS does not result from 
SN explosions in a stratified ISM, since the effective EOS 
including turbulent pressure across the local ISM models 
is close to isothermal. In the past, stiff EOSs were pro- 
posed as a way of preventing the collapse of disk galaxies 
to unphysically small radii in cosmological simulations. 
However, Li et al. (2005b) showed that such collapse was 
due to violation of the Jeans criterion (Truelove et al. 
1997) in low resolution simulations, and that fully re- 
solved, isothermal disks do not collapse. Simulations 
of full galaxies including our sub-grid model will reveal 
whether there are further effects that a stiffer equation 
of state would better reproduce. 

Limitation of the steady state assumption. — The subgrid 
model developed in this paper can be easily implemented 
in cosmological simiilations. However, our approach is 
limited by the assumption of constant SN rate; we dealt 
only with steady state systems. The model is appropriate 
if the characteristic timescale tsf over which the SN (SF) 
rate averaged over ~100 pc regions changes appreciably 
is much longer than the time it takes for the system to 
reach a steady state (^^60 Myr in our simulations). The 
assumption is valid, e.g., in merging galaxies, where tsf 
is on the order of 0.1-1 Gyr. 

7. SUMMARY 

Interstellar turbulence is thought to play a major role 
in the formation and support of molecular clouds and the 

regulation of the size, thickness, and star formation rate 
of galactic disks (Mac Low & Klcsscn 2004; Elmegrecn 
& Scalo 2004) . We examined the physical characteristics 
of the ISM by constructing high-resolution 3D models of 
a stratified ISM driven by both correlated and isolated 
SN explosions. A suite of numerical experiments, each 
with a simulated volume of (0.5 kpc)^x(10 kpc), were 
performed using an adaptive mesh code at a maximum 



spatial resolution of 1.95 pc; they included vertical grav- 
ity and parameterized heating and cooling. The simula- 
tions cover wide ranges of supernova rates and vertical 
gas column densities, with initial conditions based on the 
Schmidt-Kennicutt law (Kennicutt et al. 1998, 2007). 

We find that the turbulent velocity dispersion is in- 
versely proportional to the square root of the local den- 
sity: (Tturb oc p~^^^ on >50 pc scales. The turbulent 
pressure, perturb, in such a medium is nearly constant 
at a given scale, even though the gas density varies by 
about seven orders of magnitude. This suggests that, 
for a given SN rate, a single value of turbulent pressure 
can characterize the medium regardless of the local gas 
density. When combined with the finding in Paper I 
that >90% of the total kinetic energy is contained on 
scales below 200 pc, the two charateristics of the ISM 
are conducive to developing a subgrid model for turbu- 
lent pressure as a way to represent small-scale motions. 
One such model was developed in this paper and explic- 
itly calibrated using the local ISM simulations. 

We find that the mass-weighted velocity dispersion 
(Tturb and the simulated H I linewidth (Tv h i arc nearly 
constant across a range of 512 in SN rate (see Table 2 
for values). In other words, the appropriate equation of 
state for Pturb in gas averaged over ~100 pc regions is 
close to isothermal. In our highest supernova rate model, 
superbubble blow-outs occur, and the turbulent pressure 
on large scales is >4 times higher than the thermal pres- 
sure. 

We propose a subgrid model that naturally includes 
the effect of turbulence in large-scale simulations, at least 
as far as we understand it from our small-scale models, 
and furnishes a sound, physically motivated prescription 
for including the effect of SN feedback in cosmological 
simulations. It treats turbulence as strictly a pressure 
term, although small-scale effects on the star formation 
rate could, in principle, be included in a statistical way. 
The advantage of this type of approach is that we pre- 
scribe the subgrid model from what we know (bottom- 
up) rather than what we need in order to resolve current 
problems in cosmological simulations (top-down). 

We arc grateful to I. Goldman, Y. Li, .1. Oishi, J. Stone, 
E. Vazquez-Semadeni, and T. Thompson for stimulating 
discussions. We thank R. Cen and B. Draine for use- 
ful discussions on the scaling of the diffuse heating rate 
and J. Maron for proposing the thought experiment de- 
scribed in § 6. The software used in this work was in part 
developed by the DOE-supported ASCI/Alliance Center 
for Astrophysical Thermonuclear Flashes at the Univer- 
sity of Chicago. Computations were performed at the 
Pittsburgh Supercomputing Center and at the National 
Center for Supercomputing Applications supported by 
the NSF. 
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